Fluid Invasion in Porous Media: Viscous Gradient Percolation 
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We suggest that the dynamics of stable viscous invasion fronts in porous media depends on the 
volume capacitance of the media. At high volume capacitance, our network simulations provide 
numerical evidence of a scaling relation between the front width and its velocity. In the low volume 
capacitance regime, we derive a new effective scaling supported by network simulations and is in 
agreement with previous experiments on imbibition in paper and collections of glass beads. 

PACS numbers: 47.55.Mh, 05.40.-a, 64.60.Ak, 68.35.Fx 



Percolation theory has contributed much to our under- 
standing of immiscible fluid displacement in porous me- 
dia plEl- For slow flow, randomness in capillary forces 
dominates and invasion percolation applies. In the pres- 
ence of gravity, the local percolation probability is spa- 
tially non-uniform and gradient percolation is the stan- 
dard description Alternatively, when a viscous fluid 
displaces a non- viscous one, the pressure gradient also in- 
duces a gradient in the local percolation probability. The 
invasion front is hence expected to follow percolation ge- 
ometry as well. However, the pressure field is not known 
apriori. This problem is therefore much more challenging 
and not well-understood. 

Applying percolation theory, Xu, Yortsos and Salin 
suggested that the width of a viscous invasion front w 
depends on its propagation velocity v according to 

w~ v - K (1) 

where k ~ 0.38 in two dimensions Q. This result is 
identical to that from a Buckley-Leverett type theory of 
Wilkinson 0] although there have been other suggestions 
0, ■ Previous network simulations lead to inconclusive 
results due to the limited network sizes used 0, @ ■ In 
this paper, we report simulations at much larger scales 
using the simplest possible network models. Surprisingly, 
we observe two distinct scalings depending on the volume 
capacitance of the network defined as the volume of liquid 
that can be extracted from the porous medium locally per 
unit decrease in pressure At high capacitance, our 
simulations support the theory in Refs. 0, • At low 
capacitance, we obtain a different scaling in agreement 
with two interesting experiments on imbibition in paper 
and collections of glass beads • 
Our models are based on a simple pipe network in Ref. 
|l2| simulating the wetting of two-dimensional porous 
media with the bottom immersed into a liquid reservoir. 
We adopt a square lattice with spacing a and periodic 
boundary conditions in the lateral direction. We define a 
in-plane bond population p which is the probability that 
a bond on the square lattice is occupied by a cylindrical 
pipe of radius r = 0.25a. A fraction 1 — po of the bonds 
are left vacant. Given a pressure gradient VP along a 
pipe, the liquid flows according to Poiscuille's law for 



FIG. 1: Snapshots of fluid invasion for networks with low (a) 
and high (b) volume capacitance at in-plane pipe population 
po = 0.53. Wetted (dry) regions are shaded in dark (light) 
gray. Radii are rescaled by a factor 0.5 for clarity. Only 
connected pipes are shown. Arrows indicate instantaneous 
directions of menicus movements. Dewetting is forbidden in 
(a). In (b), one menicus recedes due to suction by other 
invaded pipes with stronger capillary forces. 



viscous flow with a flux irr VP/8/x where p. is the viscos- 
ity 0]. The atmospheric pressure is assumed to be zero 
without loss of generality. Nodes at the bottom row are 
directly connected to the liquid source and are also main- 
tained at zero pressure. The pipes are simplified models 
of complicated fluid channels. Without strictly following 
the properties of realistic pipes, we assume that the cap- 
illary pressure for an advancing meniscus inside a pipe 
follows an uniform random distribution in the interval 
[0, To] with To being a constant. The liquid pressure be- 
hind a meniscus hence varies from —To to 0. All nodes 
are volumeless and gravity is neglected. Air flow as well 
as trapping are not considered. A simple system of units 
in which a = Tq = /i = 1 is assumed. 

Two variants of the model are considered. They differ 
in the volume capacitance defined as the volume of liquid 
extracted from theporous medium locally per unit de- 
crease in pressure |9j- We first define a low capacitance 
network (Fig. da)). In this case, we forbid any reced- 
ing meniscus by assuming an ideally hysteric capillary 
pressure withstanding any tendency of dewetting. The 
liquid surface is rigid and the volume capacitance is in- 
deed zero. We also consider a high capacitance network 
(Fig. IJb)) in which receding menisci are allowed and 
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are non-hysteric. Furthermore, an extra dangling pipe 
perpendicular to the lattice is connected to every node. 
They all have length a and radius r = 0.5a. The capillary 
pressure takes a random value uniformly distributed in 
[0,0.3ro]. This enhance the volume capacitance at pres- 
sure around — 0.3r to which is the relevant range since 
it coincides with typical local pressures measured in wet- 
ted regions in our simulations. An increase (decrease) in 
the local pressure within this range activates the filling 
(draining) of many of these pipes resulting in the strong 
capacitance enhancement. We have checked that further 
increasing the capacitance by adding even more dangling 
pipes or increasing their radii give similar results. For 
simplicity, no new menisci can be generated by breaking 
any continuous column of liquid during dewetting. 

We simulate fluid invasion using Euler's method. For 
both models, the pressure field is first computed at the 
beginning of each iteration. Specifically, Poiseuille's law 
and fluid conservation at all nodes leads to Kirchoff's 
equations coupling the pressure at neighboring nodes. 
We also apply the boundary conditions that the pres- 
sure at moving and stationary menisci follows from the 
capillary pressure and the no flow condition respectively. 
The Kirchoff's equations are then solved using successive 
over-relaxation method. For the low capacitance model 
with ideally hysteric capillarity, meniscus in a partially 
filled pipe may start or stop advancement depending on 
the calculated pressure. This alters the boundary con- 
ditions and the whole calculation is repeated until the 
states of motion of all menisci are consistent with the 
pressure field. For both models, we speed up the cal- 
culation by directly computing only the pressure at the 
percolation backbone. We require that the liquid influx 
from the source agrees with the total out-flux at the mov- 
ing menisci within 1%. From the pressure, the velocities 
of all moving menisci are found. Their positions are then 
advanced over a short period in which all displacements 
are less than O.fa. 

Figure |3 plots the front width w against its velocity v 
for both models at pipe population po = 0.50, 0.53 and 
0.56. The lateral width of the lattice used is L — 1000a 
and all results have been averaged over 60 independent 
runs. Results are extracted from surface height profiles 
h(x) measured from the liquid source which mark the 
highest invaded node at lateral coordinate x. The time 
derivative of the spatial and ensemble average h gives v 
while the r.m.s. fluctuation gives w. The data fall nicely 
into straight lines in the log- log plot and verify Eq. (T{|. 
Contrary to conventional belief, there are two distinct 
sets of slopes averaging to n = 0.47(2) and 0.36(2) for the 
low and high volume capacitance networks respectively. 

To attain better insights into the problem, we now ex- 
amine closely the geometry of the invasion fronts. Fronts 
at various locations and widths will be compared. We 
hence define a rescaled distance z from the center of 
a front by z = (y — h)/w where y is the vertical co- 
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FIG. 2: Log-log plot of the invasion front width w against 
velocity v. The slopes of the dashed lines are -0.47 and - 
0.36 respectively. Inset: log-log plot of w against vlnw. The 
dashed line has a slope -0.53. 



ordinate. An important quantity is the local satura- 
tion S(y) defined as the fraction of in-plane pipes which 
are completely invaded at height y. Figure 03a) plots 
S(y)/w D ~ 2 against z using the same data as in Fig. Q] 
where D = 91/48 is the fractal dimension 0. Values 
recorded near the beginning and the end of each set of 
simulations are plotted. Specifically, time t goes from 
1.8 x 10 5 to 2.5 x 10 8 corresponding to w in between 3.1 
and 32. All 12 curves collapse nicely into a single curve 
implying the scaling form 



S(y) 



n-2 



fs(z) 



(2) 



where fs is a rescaled function. This verifies the perco- 
lation geometry of the fronts with w being the correla- 
tion length 0. Similarly, we define the local percolation 
probability p(y) as the probability that an in-plane pipe 
connected to an invaded node is completely filled. Fig- 
ure|3fb) plots (p(y) —p c )/w~ 1 / L ' against z where v = 4/3 
and p c = 1/2 are respectively the correlation length expo- 
nent and the percolation threshold. Good data collapse 
is again observed supporting 



p(y) ~Pc = w 1,w f p {z). 



(3) 



Data collapses for pressure in Fig. I^c) are most in- 
teresting and will be explained after examining the re- 
lation between the local pressure and the local percola- 
tion probability. During invasion, the percolation prob- 
ability at any point increases as pipes with ever lower 
capillary pressure are filled. Given a local percolation 
probability p, pipes with capillary pressure as low as 
F(p) = (1 — p/po)T Q are expected to be filled. There 
is hence at least one instance when the pressure climbs 



3 



0.7 
0.6 
0.5 

Q 0.4 

3 

w 0.3 

0.2 

0.1 


0.3 
0.2 

> 0.1 

? 3 ° 
^ -0.1 

s£ -0.2 
a. 

-0.3 
-0.4 
-0.5 



0.2 


-0.2 
-0.4 
-0.6 
-0.8 



Low vol. cap. 

□ p = 0.50 
a p = 0.53 
© p = 0.56 

High vol. cap. 

• p = 0.50 
» p = 0.53 

♦ p = 0.56 



< 




FIG. 3: Rescaled plots of local saturation S (a), percolation 
probability p (b), and pressure P (c) against rescaled coordi- 
nate z = (y — h)/w. Note that P requires different rescaling 
forms for the low and high volume capacitance cases. 



up to instantaneous maximum value —Tip). Figure 01 
compares the local pressure P with —Tip). Only results 
for p$ — 0.53 are shown but other values give similar 
findings. 

At high volume capacitance, Fig. 0] shows that the 
pressure simply coincides with its upper bound, i.e. 



(4) 



(5) 



P = -T{p) 
Letting P c = — T(p c ), it implies 

p (y) - p c ~ p(y) - Pc 



which has been widely used in the literature [_ 
and is valid for gradient percolation |j| . However, for vis- 
cous gradient percolation, Eq. (@J is indeed non-trivial 
and results from the strong damping of pressure fluc- 
tuations. Local pressure depends on the capillary pres- 
sure of the pipes under invasion. Pipes with capillary 
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FIG. 4: Plot of the average pressure P against the instanta- 
neous maximum pressure —Tip) expected from the percola- 
tion probability p. 



pressure in the range [r(p),rn] are to be invaded. The 
pressure is hence dominated by those with the small- 
est capillary pressure Tip) which takes the longest time 
to wet. Moreover, we observe that when a pipe with a 
higher capillary pressure is invaded, it quickly draws liq- 
uid from neighboring pipes, especially perpendicular ones 
with capillary pressure close to Tip). There is hence only 
a brief and highly localized pressure disturbance which 
has little impact on the average pressure. After com- 
pleting the invasion, these neighbors are slowly refilled 
directly from the liquid reservoir and maintain the pres- 
sure roughly at —T(p) again. This pressure stabiliza- 
tion mechanism leads to Eq. J3J|. Then, Eqs. J3J and 
© imply P(y) - P c = w~ 1/u 'fp(z). We therefore plot 
(P(y) — P c )/w~ 1 l v against z for the high capacitance 
case and is shown in Fig. Etc) using filled symbols. We 
have taken v' = 1.5 for the best data collapse which is in 
reasonable agreement with the exact value 4/3. The dis- 
crepancy inherits from the slight deviation from Eq. Q 
as observed in Fig. 0]and results from imcomplete damp- 
ing of the pressure fluctuations due to the finite volume 
capacitance of the network. The pressure difference A.Ph 
across the front is hence 



AP 



H 



(6) 



At low capacitance, Fig. 0] implies P < —Tip) in- 
stead. Pipes with capillary pressure larger than T(p) can 
now significantly pull down the pressure when invaded. 
A new analytic description of the pressure field is now 
presented. We first quantitatively define the front region 
by z a < z < Zb- We take z a — — 2 since it marks the 
bottom of the lattice for po — 0.50 and = 2 as it 
is where the saturation practically vanishes. We denote 
quantities evaluated at the boundaries by the subscripts 
a and b. For z < z a corresponding to the bulk region, 
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good connectivity suppresses pressure fluctuations and 
Eq. holds approximately. In particular, the pressure 
at z a is P a = —T(p a ). In contrast, the pressure fluc- 
tuates strongly at Zf, between — r and — T(jpb) as pipes 
with capillary pressure in the range [r(p ),Fo] are in- 
vaded. With dewetting in other pipes forbidden, they 
draw liquid directly from the water reservoir. This leads 
to strong delocalized pressure fluctuations. The inva- 
sion of a pipe with capillary pressure P cap takes a period 
t ~ {P C ap + -Pa) -1 - A characteristic pressure P at z 
hence follows from the time averaged capillary pressure. 
Straight-forward algebra gives P — P a — APl where 



AP L = [T -r(p 6 )]/ln 



P, 



T(p b ) + P a 



(7) 



Here APl represents a characteristic pressure difference 
across the front. We thus suggest P(y)—P a = APl/p(z). 
It is verified by the reasonable data collapse in Fig. |3{c) 
which plots (P(y) — P q )/APl against z for the high 
capacitance case using open symbols. We have used 
Pa,b — Pc + w^ 1 /" fp(z a ^) from Eq. (JSJ) and further as- 
sumed f p {z a ) = 0.15 and f p (zb) = —0.5 which are con- 
sistent with values read directly from Fig. &[b). 

Now, we can calculate the exponent k defined in Eq. 
JIJ. The permeability of the front is k ~ u; - ^/" -1 where 
H is the conductivity exponent 0] • The liquid flux across 
the front is then kAPji,L for high and low capacitance 
respectively. However, it should equal vS as well. Using 
also Eq. (S), we obtain v ~ wV/ v + D - x AP l ,h- For the 
high capacitance model obeying Eq. © , we arrive at Eq. 
JTJ) with 

k = v/(1 + [i + v(D-1))~0.38 (8) 

derived previously in Refs. 0,0,0- Our network simula- 
tions in the high capacitance regime lead to n ~ 0.36(2) 
giving a direct numerical support to Eq. ijSJl. 

For the low capacitance case, we apply Eq. (7J| which 
simplifies at large w to APl ~ (lniu) -1 . Therefore, Eq. 
||TJ is replaced by 



u,-V« 
law 



with a new exponent 



k = v/(ji + v(D - 1)) -0.53 



(9) 



(10) 



The scaling thus admits a logarithmic correction. This is 
verified by the log- log plot of w against v In w in the inset 
of Fig. which gives straight lines with the expected 
slope -0.53(2) for w > 8. A naive measurement of k from 
a log-log plot of w against v disregarding the correction 
should lead to an effective exponent K e in between 0.38 
and 0.53. Using Eq. JJJ, we obtain n e ~ 0.46 in good 
agreement with the value 0.47(2) directly measured from 
our simulations and 0.48 from imbibition experiments in 
Refs. [ill ED- 



We have considered extreme values of the volume ca- 
pacitance. In general, networks with significant pressure 
fluctuations extending beyond or well- within the invasion 
front are at the low or high capacitance limit respectively. 
Networks with intermediate capacitance are therefore ex- 
pected to crossover to the high capacitance regime at the 
large front width limit. However, from experiments, pres- 
sure fluctuations in the form of Haines jumps propagate 
far beyond the invasion front 9J. It is also known that 
wetting and drainage in porous media are highly hysteric. 
These further establish that the new low-capacitance con- 
dition is the correct experimentally relevant re gim e. In 
addition, similar to the original model in Ref. [l2j, our 
low capacitance model can also reproduce the rich dy- 
namical features of the invasion front in Ref. ^3 with 
reasonable accuracy at po = 0.53 and will be reported 
elsewhere. 



In conclusion, we have derived a new effective scaling 
relation based on viscous gradient percolation for fluid 
invasion fronts in porous media with low volume capac- 
itance. This regime is characterized by realistic features 
including highly hysteric capillary flow and long-range 
propagation of pressure disturbances. The scaling is sup- 
ported by large scale network simulations and agrees with 
previous experiments. We also show numerically that a 
well-known scaling theory applies only at high volume 
capacitance which is not supported experimentally. 

We thank V.K. Horvath, L.M. Sander, J.Q. You and 
F.G. Shin for helpful discussions. This work was sup- 
ported by HK RGC Grant No. PolyU-5191/99P. 
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